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TECHNICAL PAPER 


APPLICATIONS OF FEM AND BEM IN TWO-DIMENSIONAL 
FRACTURE MECHANICS PROBLEMS 

I. INTRODUCTION 


The finite element method (FEM) and boundary element method (BEM) are widely employed in 
fracture mechanics. 1-6 Both methods are employed for fracture mechanics problems, but neither method 
is effective and accurate for all fracture problems. Judging the effectiveness and accuracy of FEM or 
BEM requires practical experience in the field. Traditionally, the stress analysis has used the FEM. 
However, the fracture mechanics problems have crack edge boundary and singular points at the crack tip 
which are inconvenient to mesh because full domain discretization is required. In addition, around crack 
tips, extremely fine subdivisions are required to achieve reasonable accuracy. It takes quite a length of 
time for an engineer to build and check such a model. It has been proposed that these disadvantages are 
overcome by the BEM since the elements are placed only on the boundary of the domain. However, this 
method also has difficulty in mathematical degeneracy and the treatment of nonlinearities. A number of 
different approaches have been developed to address the problem of degeneracies. Also, in order to 
overcome the disadvantages of these methods, coupled finite element — boundary element methods have 
been developed. 7 

The objective of this paper is not to make a judgment on these methods but, to show how the 
FEM and BEM can be employed in fracture mechanics problems and also to find how both methods can 
be used in practical applications, such as analyzing component failure in aerospace structures caused by 
fracture and fatigue. 


H. FINITE ELEMENT METHOD 


Since several excellent and comprehensive works explaining basic fracture mechanics theory are 
readily available in the literature, 7-12 a comprehensive review is not made here. However, selected fea- 
tures of the fundamental theory are reviewed. A critical step in modeling fracture problems is the selec- 
tion of the type of element used to represent the singular region surrounding the crack tip. If regular (i.e., 
nonsingular) elements with polynomial displacement functions are selected, the convergence rate, based 
on element size, is very slow. 13 Moreover, the convergence rate at a singularity cannot be increased by 
merely selecting higher order polynomial elements. The ANSYS six-node triangular element with mid- 
side nodes at the quarter points was selected for the crack tip since it accurately models the singularity. 14 
The rest of the domain was meshed with standard six-node triangular elements. The development of the 
singular isoparametric element is well documented in the literature. 1 6 15-19 The stress intensity factors 
were calculated from the nodal displacements, using the same techniques as the BEM models discussed 
in section III. The crack-tip near-field displacements and stresses for mode I behavior (the opening 
mode) and mode II (the shearing mode) are expressed in equations (1) and (2) in the crack-tip coor- 
dinates which are shown in figure l. 10 
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in which K\ and Kjj are the stress intensity factors associated with mode I and mode II, respectively, of 
linear elastic fracture mechanics (LEFM). cr x , o y , and t 'xy are the stress components in the rectangular 
(x,y,z) coordinate system; u and v are the displacements; r and 0are local polar coordinates; G is the 
shear modulus; and vis Poisson’s ratio. The variable jcis a conversion factor between plane strain and 
plane stress, where 


k = 3-4 v for plane strain 

3-v ^ 

= y-j-^ for plane s tres s . 

The evaluating procedure for Kj and Kji from the finite element nodal displacements is not new 17 20 21 23 
and neglecting higher order terms. 



where only the absolute value of displacement is considered. 


Thus, the full opening and shearing displacements across the faces of the crack are obtained as 
follows: 



Av = — 
2G 



( 5 ) 


2 



The method to extract AT/ from mode I displacements is described here since the procedure is 
identical for AT//. In this case, the behavior of either v(r) in equation (4) or Av(r) in equation (5) is 
approximated. On the crack face, it can be written 


£ or & = A+Br , (6) 

Vr Vr 

where v and Av are the approximations to the crack-face displacements, and A and B are constants 
determined from a linear curve fit of nodal displacements. All displacements are relative to the crack-tip 
node. Substituting the values of v or Av and r for nodal points along a radial line emanating from the 
crack tip allows a plot of equation (6) against radial distance r to be drawn as shown in figure 2. Then by 
discarding the results for points very close to the crack tip, the solutions can be extrapolated to r = 0. 
Once A and B are determined, then the limit is found as r-» 0. 


lim 
r-> 0 


x 

Vr 


or 


lim Av _ a 


(7) 


Solving first equation (4) and then equation (5) for Kj and combining with equation (7), 

K[ = 2(7 for half-crack models 
( 1 +k) 

K[ = GJ2kA for full-crack models . 

(1 + k) 


( 8 ) 


Using these steps, a single-edge crack and an inclined-edge crack in an elastic plate were solved and the 
results are described in the following subsections. 

A. Elastic Analysis of a Single-Edge Crack 

The two-dimensional plate geometry is shown in figure 3, and the plane strain assumptions are 
used in this example. The model is symmetric about the y = 0 axis, and the x-displacement at x = 0 and 
y = 0 was restrained to prevent rigid body translation. The dimensions used were alb = 1 and da = 6, 
where a is the crack length, and b and c are plate dimensions. A uniform tensile stress, cr, is applied at 
the ends of the plate. Material properties are Young’s modulus of 10xl0 6 lb/in 2 and Poisson’s ratio of 
0.333. The theoretical stress intensity factor, AT/, associated with the crack-opening mode is expressed as 

Kj = afna-fialb) , (9) 

where f(a/b) is a geometric correction factor. For a/b, a nearly exact solution is f(a/b ) = 2.82. 10 It is 
convenient to rewrite equation (9) in dimensionless form: 


Ki = 


K, 

G-fita 


=f(a/b) . 


( 10 ) 


This leads to AT/ = 2.82 as the exact solution to compare to the finite element results. In the model, six- 
node triangular elements are used except for eight singular elements placed at the tip of the crack. Also, 
only half the plate is modeled because of symmetry along the direction of the crack. The analysis results 
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are given in figure 4 and table 1 and indicate favorable agreement with the exact solution. Also, since 
most real engineering fracture problems involve more than a single mode of fracture, the next problem 
attempted was an inclined-edge cracked plate. 

B. Elastic Analysis of an Inclined Crack 

An elastic two-dimensional plate, weakened by the presence of an inclined edge crack as shown 
in figure 5, was selected to extract the Kj and Ku stress intensity factors under plane strain assumptions. 
The boundary conditions are v = 0 along the x-axis and u = 0 at the point (0,0). The dimensions used 
were aAv = 0.3 and L/w = 5.0 where a is crack length, w is plate width, and L is plate length. The crack is 

inclined 45°. The plate is subjected to a tensile stress (o) at its ends that will induce crack face opening 
(mode I) and crack face shearing (mode II). Since no planes of symmetry exist, a full model was gener- 
ated with six-node triangular elements and singular elements placed in the crack-tip region. The analysis 
results are shown in figure 6 and table 1. Having the crack faces defined as shown in figure 2, the values 

of Ki = 1.56 and Ku = 0.786 were obtained. The theoretical values reported by Wilson 22 are 1.58 and 
0.782, respectively. 


HI. BOUNDARY ELEMENT METHOD 


The BEM has been developed by a number of researchers, and recently a great deal of interest 
has been focused on this method as a complement to the FEM. 24 Although this technique offers accurate 
results with significant savings in time and/or cost, the range of practical applications is currently not as 
wide as the finite element technique. One area that has been identified as being well-suited for boundaty 
element analysis is fracture mechanics since it only requires discretization of the surface rather than the 
volume. The simplified theory of the BEM follows: 25 

Take a component with volume £2 and boundary T (fig. 7). The component is subjected sepa- 
rately to two load cases: 

Load case 1: Loads t, Displacements u 
Load case 2: Loads t*. Displacements w*. 

The reciprocal theorem states that the work done by the loads from load case 1 on the displace- 
ments from load case 2 is equal to the work done by the load case 2 on the displacements from 
load case 1. To find the work, we simply integrate over the body, so we can write 

= \j *' u • ( 11 ) 

Assuming that there are no body loads, the loading is confined to the boundary T, so the volume 
integrals may be rewritten as surface integrals. Body forces may still be considered in the BEM 
by using a Green’s function to transform the volume integral term to a boundary integral. 
Rewriting as surface integrals. 


j r t-u*dT = ^t*-udT . 


( 12 ) 
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Now we specify that load case 1 represents the actual loading for which the analysis is per- 
formed. Load case 2, still being arbitrary, can be chosen to be anything we want to make the 
equations easier to handle. In practice, it turns out to be best to apply a point force at the node 
locations (like finite elements, boundary elements have nodes at which the equations are 
formed). We know the terms t* and u* from classical solutions. It is usual to call u* “essential” 
and t* “natural” boundary conditions. 

Integrating the terms t and u numerically over the elements, a system of equations is formed in 
terms of the nodal displacements and loads 

[H] [U] = [G] [7] , (13) 

where [if] and [G] are two NxN matrices and [U] and [7] are vectors of length N. 

The square matrix [77] is filled with the integrals of the t* terms, and the square matrix [G] is 
filled with the integrals of the u* terms. The vector [71 contains the unknown loads and the vec- 
tor [U] contains the unknown displacements. This leaves us with TV equations, but 2N unknowns. 

However, since displacements (i.e., u*) on Ti and tractions (i.e., t*) on T 2 are known, there are 
only N unknowns in the system of equation (13). To introduce these boundary conditions into 
equation (13), one has to rearrange the system by moving columns [77] and [G] from one side to 
the other. Once all unknowns are passed to the left-hand side, one can write, 

[A][X\ = [F] , (14) 

where [X] is a vector of unknowns u and t boundary values. [F] is now all known. 

This gives a familiar looking equation which may be solved by any standard matrix equation 
solver. Once the boundary solution has been found, it is possible to calculate any internal value 
of u or its derivatives. The values of u are calculated at any internal point “j” using the formula 
which can be written as. 


n‘ = J r t • u*dT — j u • t*dT . (15) 

Notice that now the fundamental solution is considered to be acting on an internal point “i” and 
that all values of u and t are already known. The process is then one of integration (usually 
numerically). The stress and displacement at any given point internal to the material can be found 
by a similar integration procedure. 

Although the theory underlying the BEM is relatively complicated, model creation is relatively 
simple because only the boundary needs to be defined. The details of modeling crack tips with boundary 
elements varies according to the software chosen. The Boundary Element Analysis System (BEASY) 
was used for this analysis. The following description is based on the BEASY code. In BEASY models, 
line elements are used for two-dimensional problems. The points which define the geometry of these ele- 
ments are called “mesh points.” Unlike the FEM, these mesh points do not necessarily coincide with the 
“nodes,” at which the equations are formed for the results. It is necessary to make the distinction 
between mesh points and nodes. Mesh points define only the geometry, and the nodes define the 
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location where the values of displacement, temperature, etc., on the element are calculated. In most 
instances, the nodes are located at the mesh points for continuous elements. But when discontinuous 
elements are needed (as is the case for a discontinuous stress field at a crack tip), the nodes are located 
away from the mesh point. The number of nodes on an element varies according to the user’s 
requirements. Quadratic elements were used in this study. The nodal results were calculated and then 
extrapolated to the mesh points. Mesh points are always at the middle and ends of two-dimensional 
elements and are shared wherever possible between adjacent elements. The extrapolation of results to a 
mesh point done by fitting a quadratic curve through the three nodal results can, therefore, be done from 
two or more elements. These mesh point results are, therefore, the mean values of the results obtained 
from the different elements connected to a mesh point. 

As shown in figures 8 and 9, the models are meshed with graded element sizes. No special ele- 
ments need to be defined at a crack-tip, but smaller elements are needed. In the region approaching a 
crack tip, the stress is varying very rapidly, and it is, therefore, important to use small elements. Use of 
such small elements throughout the model, though, would create too many elements to be run prac- 
tically, so they should be graded out to larger elements in other parts of the model. However, if run time 
is not a dominating factor, then users can define very small elements at the crack tip. It should also be 
noted that a boundary element model can only model accurately if the elements are not required to 
model stress variation of a higher order than the elements themselves. For example, a quadratic element 
will be accurate only if the stress does not vary as a cubic function (or higher) over the length of the 
element. 

The stress intensity factors are calculated using Irwin’s equations 26 from both stress and dis- 
placement. The latter have been found to give consistently more accurate results. This is perhaps due to 
the fact that stress becomes singular at the crack-tip, whereas displacement does not. Notice that these 
equations are not valid at the crack tip. By extrapolating these values to the crack tip itself, the user can 
obtain an accurate estimate of the true stress intensity factor. In this study, the stress intensity factors 
were normalized as defined in section II. The same models solved in section II were modeled by 
BEASY and the results were obtained as follows. 

A. Elastic Analysis of a Single-Edge Crack 

The exact same model used for the finite element analysis shown in figure 3 was used for the 
boundary element analysis. The model contained 65 line elements. Only the top half of the plate was 
modeled due to the symmetry. Figure 8 illustrates the BEASY model of this problem. Following the 
analysis, a plot of the deformed shape was made, and the crack opening was observed as shown in this 

figure. The normalized Kj factor was obtained by the method described above, and the value of A) was 
2.88. The result was then compared to the solution given in reference 10. The difference was 
approximately 2.1 percent. The theoretical value reported by Paris and Sih 10 is 2.82. 

B. Elastic Analysis of an Inclined Crack 

A crack developing from the edge of a plate at 45° to the plate edge is shown in figure 5. The 
whole plate was modeled using BEM because of no symmetry. In this case, the shear mode (mode II) 
acts in addition to the crack opening mode (mode I), therefore both Kj and Kjj were calculated. The 

model consisted of 1 14 elements, and gave the values of Kj = 1.57 and Kjj = 0.794. The results were then 
compared to those obtained by Wilson. 22 The difference was less than 2 percent as shown in table 1. The 
theoretical values reported by Wilson 22 are 1.58 and 0.782, respectively. The deformed shape on the 
original shape is shown in figure 9. 
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IV. RESULTS AND CONCLUSIONS 


Application of the FEM and BEM is just beginning for fracture mechanics analysis in the 
aerospace industry. Traditional methods/handbook solutions are still the predominate techniques used. 
However, these traditional techniques have limitations to simulate the practical cases. Therefore, the 
FEM and BEM have been considered and developed to resolve the limitation of traditional methods. 
However, the FEM is expensive in both manpower and computer costs, and cost saving alternatives are 
always being sought. The BEM has been considered to be one of them. 

The single edge cracked plate and inclined edge cracked plate were analyzed with both the 
BEASY boundary element code and the ANSYS finite element code to see the relative strengths and 
weaknesses of both methods. Stress intensity factor calculations were computed on all the finite element 
models and boundary element models using the displacement extrapolation method. Final calculations 
for the stress intensity factors for the r fit are given in table 1. The boundary element K predictions were 
within 2 percent of the finite element K predictions. The agreement between the boundary element 
model and the finite element model is encouraging considering the difference in number of degrees of 
freedom (DOF) (1,250 with FEM to 266 with BEM for single edge crack case, 2,226 with FEM to 536 
with BEM for inclined edge crack case). However, from a simple comparison of the solution times, the 
BEM (BEASY code) does not show more efficiency than the FEM (ANSYS code) as can be seen in 
table 1. For this comparison, a VAX/785 computer system was used. The most probable reason for this 
is that the boundary element model uses a “full” matrix without the banded symmetry common to the 
FEM. In other words, the relationship between matrix size, fullness, and the central processing unit 
(CPU) time to invert and solve the equations might be the reason for the efficiency in computing time. 
The FEM was more accurate than the BEM for the single edge crack problem. However, the Kj from the 
BEM was more accurate than that from the FEM, and the Kji from the FEM was more accurate than that 
from the BEM for the inclined edge crack problem. The most significant observation was that both 
methods provided a good correlation with the theoretical predictions, but about the same amount of 
modeling time was required to achieve favorable agreement with the exact solutions using both methods. 
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Table 1. Comparison of FEM and BEM. 


Method 

Case 

Single Edge 

Inclined Edge 

FEM 

Number of Elements 

292 

526 

Sol. Time CPU 

190.15 s 

384.85 s 

K, 

2.83 

1.565 

Percent A 

0.4 percent 

-0.95 percent 

Kn 

— 

0.7857 

Percent A 

— 

0.47 percent 

DOF 

1,250 

2,226 

BEM 

Number of Elements 

65 

114 

Sol. Time CPU 

254.30 s 

727.50 s 

K, 

2.88 

1.570 

Percent A 

2.1 percent 

-0.6 percent 

K n 

— 

0.794 

Percent A 

— 

1.5 percent 

DOF 

266 

536 

Ref. 

K, 

2.82 [10] 

1.58 [22] 

Kn 

— 

0.782 [22] 



Figure 1. Local coordinates measured from a two-dimensional crack tip. 
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Deformed Shape 


Original Shape 


X 



Figure 9. Deformed shape on original shape. 
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